Driving Sandpiles to Criticality and Beyond 



O 
(N 

o 

CD 



O 
> 

o 

m 

O 



X 



I Anne Fey} ^ ILionel Levinep and IDavid B. Wilsonf 

'Delft Institute of Applied Mathematics, Delft University of Technology, The Netherlands 
"^Department of Mathematics, Massachusetts Institute of Technology, Cambridge, MA 02139, USA 
■^Microsoft Research, Redmond, WA 98052, USA 
(Dated: December 16, 2009; revised March 4, 2010) 

A popular theory of self-organized criticality relates driven dissipative systems to systems with 
conservation. This theory predicts that the stationary density of the abehan sandpile model equals 
the threshold density of the fixed-energy sandpile. We refute this prediction for a wide variety of 
underlying graphs, including the square grid. Driven dissipative sandpiles continue to evolve even 
after reaching criticality. This result casts doubt on the validity of using fixed-energy sandpiles to 
explore the critical behavior of the abelian sandpile model at stationarity. 

PACS numbers: 64.60.av, 45.70.Cc 



A subset of the vertices are distinguished as sinks: they 
absorb particles but never topple. A single time step con- 
sists of adding one particle at a random site, and then 
performing topplings until each non-sink vertex has fewer 
particles than its degree. The order of topplings does not 
affect the outcome [l^j. The set of topplings caused by 
addition of a particle is called an avalanche. 

Avalanches can be decomposed into a sequence of 
"waves" so that each site topples at most once during 
each wave. Over time, sandpiles evolve toward a station- 
ax^ state in which the waves exhibit power-law statistics 
13] (though the full avalanches seem to exhibit multifrac- 
tal behavior IJ, ll5|). Power-law behavior is a hallmark 
of criticality, and since the stationary state is reached 
apparently without tuning of a parameter, the model is 
said to be self-organized critical. 

To explain how the sandpile model self-organizes to 
reach the critical state, Dickman et al. [H, 0] introduced 
an argument which soon became widely accepted: see, for 
example, [l^ Ch. 15.4.5] and [l7l - [l9| . Despite the appar- 
ent lack of a free parameter, they argued, the dynamics 
implicitly involve the tuning of a parameter to a value 
where a phase transition takes place. The phase tran- 
sition is between an active state, where topplings take 
place, and a quiescent "absorbing" state. The param- 
eter is the density, the average number of particles per 
site. When the system is quiescent, addition of new par- 
ticles increases the density. When the system is active, 
particles are lost to the sinks via toppling, decreasing 
the density. The dynamical rule "add a particle when 
all activity has died out" ensures that these two density 
changing mechanisms balance one another out, driving 
the system to the threshold of instability. 

To explore this idea, DMVZ introduced the fixed- 
energy sandpile model (FES), which involves an explicit 
free parameter ^, the density of particles. On a graph 
with TV vertices, the system starts with (N particles at 
vertices chosen independently and uniformly at random. 
Unlike the driven dissipative sandpile described above, 
there are no sinks and no addition of particles. Subse- 



In a widely cited series of papers [1H5|j Dickman, 
Mufioz, Vespignani, and Zapperi (DMVZ) developed a 
theory of self-organized criticality as a relationship be- 
tween driven dissipative systems and systems with con- 
servation. This theory predicts a specific relationship 
between the abelian sandpile model of Bak, Tang, and 
Wiesenfeld Q, a driven system in which particles added 
at random dissipate across the boundary, and the cor- 
responding "fixed-energy sandpile," a closed system in 
which the total number of particles is conserved. 

After defining these two models and explaining the 
conjectured relationship between them in the DMVZ 
paradigm of self-organized criticality, we present data 
from large-scale simulations which strongly indicate that 
this conjecture is false on the two-dimensional square lat- 
tice. We then examine the conjecture on some simpler 
families of graphs in which we can provably refute it. 

Early experiments 0] already identified a discrepancy, 
at least in dimensions 4 and higher, but later work fo- 
cused on dimension 2 and missed this discrepancy (it is 
very small). Some recent papers (e.g., 0]) restrict their 
study to stochastic sandpiles because deterministic sand- 
piles belong to a different universality class, but there 
remains a widespread belief in the DMVZ paradigm for 
both deterministic and stochastic sandpiles [^, [13] ■ 

Despite our contrary findings, we believe that the cen- 
tral idea of the DMVZ paradigm is a good one: the dy- 
namics of a driven dissipative system should in some way 
reflect the dynamics of the corresponding conservative 
system. Our results point to a somewhat different re- 
lationship than that posited in the DMVZ series of pa- 
pers: the driven dissipative model exhibits a second-order 
phase transition at the threshold density of the conser- 
vative model. 

Bak, Tang, and Wiesenfeld @ introduced the abelian 
sandpile as a model of self-organized criticality; for math- 
ematical background, see [11[. The model begins with a 
collection of particles on the vertices of a finite graph. 
A vertex having at least as many particles as its degree 
topples by sending one particle along each incident edge. 
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quently the system evolves through topphng of unstable 
sites. Usually the parallel toppling order is chosen: at 
each time step, all unstable sites topple simultaneously. 
Toppling may persist forever, or it may stop after some 
finite time. In the latter case, we say that the system 
stabilizes; in the terminology of DMVZ, it reaches an 
"absorbing state." 

A common choice of underlying graph for FES is the 
n X n square grid with periodic boundary conditions. It 
is believed, and supported by simulations [l^ , that there 
is a threshold density (^c, such that for C < Cc, the system 
stabilizes with probability tending to 1 as n — cxd; and 
for C > with probability tending to 1 the system does 
not stabilize. 



THE DENSITY CONJECTURE 

For the driven dissipative sandpile on the n x n square 
grid with sinks at the boundary, as n oo the stationary 
measure has an infinite- volume limit |21J , which is a mea- 
sure on sandpiles on the infinite grid Z^. One gets the 
same limiting measure whether the grid has periodic or 
open boundary conditions, and whether there is one sink 
vertex or the whole boundary serves as a sink (21| (see 
also [22I for the corresponding result on random span- 
ning trees). The statistical prope rties o f this limiting 
measure have been much studied 
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Grassberger 

conjectured that the expected number of particles at a 
fixed site is 17/8, and it is now known to be 17/8± 10^^^ 
[2^. We call this value the stationary density of Z^. 

DMVZ believed that the combination of driving and 
dissipation in the classical abelian sandpile model should 
push it toward the threshold density Cc of the fixed- 
energy sandpile. This leads to a specific testable pre- 
diction, which we call the Density Conjecture. 

Density Conjecture On the square grid, — 

17/8. More generally, (^c — Cs- 

Vespignani et al. ,4] write of FES on the square grid, 
"the system turns out to be critical only for a partic- 
ular value of the energy density equal to that of the 
stationary, slowly driven sandpile." They add that the 
threshold density Cc of the fixed energy sandpile is "the 
only possible stationary value for the energy density" of 
the driven dissipative model. In simulations they find 
Cc = 2.1250(5), adding in a footnote "It is likely that, 
in fact, 17/8 is the exact result." Other simulations to 
estimate Cc also found the value very close to 17/8 [HQ- 

Our goal in the present paper is to demonstrate that 
the density conjecture is more problematic than it first 
appears. Table U presents data from large-scale simula- 
tions indicating that Cc(^^) is 2.125288 to six decimals; 
close to but not exactly equal to 17/8. 

In each trial, we added particles one at a time at uni- 
formly random sites of the n x n torus. After each ad- 



n 


trials 


estimate of ^c(2^^^) 


64 


22s 


2.1249561 ± 0.0000004 


128 


226 


2.1251851 ± 0.0000004 


256 


224 


2.1252572 ± 0.0000004 


512 


222 


2.1252786 ± 0.0000004 


1024 


220 


2.1252853 ± 0.0000004 


2048 


2IS 


2.1252876 ± 0.0000004 


4096 


2I6 


2.1252877 ± 0.0000004 


8192 


2l4 


2.1252880 ± 0.0000004 


16384 


2^^ 


2.1252877 ± 0.0000004 





2.125288 




/ 

/ 










; 


2.125000000000 




n 

H \ \ \ \ \ \ \ h> 



64 128 256 512 1024 2048 40% 8192 16384 



TABLE I: Fixed-energy sandpile simulations on n x n tori Z^j . 
The third column gives our empirical estimate of the thresh- 
old density Cc(Z^) for Z^. The standard deviation in each 
of our estimates of Cc(^n 

) is 4 X 10^^. To six decimals, 
the values of Cc(^2048)) • ■ • i CcC^iessi) a-re all the same. The 
data from n — 64 to n = 16384 are well approximated by 
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2.1252881 ± 3 X 10" 



(0.390 ± 0.001)n- 



shown in the graph. (The error bars are too small to be visi- 
ble, so the data are shown as points.) The rapid convergence 
is due in part to periodic boundary conditions. We conclude 
that the asymptotic threshold density Cc(Z^) is 2.125288 to 
six decimals. In contrast, the stationary density Cs(^^) is 
2.125000000000 to twelve decimals. 



dition, we performed topplings until either all sites were 
stable, or every site toppled at least once. For determin- 
istic sandpiles on a connected graph, if every site topples 
at least once, the system will never stabilize 2^ 2^ . We 
recorded m/n^ as an empirical estimate of the threshold 



density Cc(Z^), where m was the maximum number of 
particles for which the system stabilized. We averaged 
these empirical estimates over many independent trials. 

We used a random number generator based on the Ad- 
vanced Encryption Standard (AES-256), which has been 
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found to exhibit excellent statistical properties 
Our simulations were conducted on a High Performance 
Computing (HFC) cluster of computers. 



PHASE TRANSITION AT Cc 

We consider the density conjecture on several other 
families of graphs, including some for which we can de- 
termine the exact values Cc and Cs analytically. 

Dhar ^1^ defined recurrent sandpile configurations and 
showed that they form an abelian group. A consequence 
of his result is that the stationary measure for the driven 
dissipative sandpile on a finite graph G with sinks is the 
uniform measure on recurrent configurations. The sta- 
tionary density Cs{G) is the expected total number of 
particles in a uniform random recurrent configuration, 
divided by the number of non-sink vertices in G. 

The threshold density Cc and stationary density Cs for 
different graphs is summarized in Table |lTl The only 
graph on which the two densities are known to be equal 
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graph 







Z 


1 


1 


I? 


17/8 = 2.125 


2.125288... 


bracelet 


5/2 = 2.5 


2.496608 . . . 


flower graph 


5/3 = 1.666667 . . . 


1.668898 . . . 


ladder graph 


1 - ^ = 1.605662.. 


1.6082 . . . 


complete graph 


1/2 X n + 0{^) 


1 X n — 0{^n log n) 


3-regular tree 


3/2 


1.50000. . . 


4-regular tree 


2 


2.00041. . . 


5-regular tree 


5/2 


2.51167. . . 



TABLE II: Stationary and threshold densities for difl^erent 
graphs. Exact values are in bold. 



is Z [17|, |l8j, |27[ . On all other graphs we examined, with 
the possible exception of the 3-regular Cayley tree, it 
appears that Cc 7^ Cs- 

Each row of Table [TT] represents an infinite family of 
graphs Gn indexed by an integer n > 1. For example, 
for we take Gn to be the nxn square grid, and for the 
regular trees we take G„ to be a finite tree of depth n. As 
sinks in G„ we take the set of boundary sites Gn \ G„_i 
(note that on trees this corresponds to wired boundary 
conditions). The value of Cs reported is lim„_j.oo Cs(G„). 

The exact values of for regular trees (Bethe lattices) 
were calculated by Dhar and Majumdar [sij. The cor- 
responding values of we report come from simulations 
[32| . We derive or simulate the values of Cs and for the 
bracelet, flower, ladder, and complete graphs in [32 1. 

As an example, consider the bracelet graph _B„, which 
is a cycle of n vertices, with each edge doubled (see Fig- 
ure [J). A site topples by sending out 4 particles: 2 to 
each of its two neighbors. One site serves as the sink. To 
compare the densities and Cs, we consider the driven 
dissipative sandpile before it reaches stationarity, by run- 
ning it for time A. More precisely, we place An particles 
uniformly at random, stabilize the resulting sandpile, and 
let Pn(A) denote the expected density of the resulting sta- 
ble configuration. In the long version of this paper 32 1 
we prove 



Theorem 1 (;32]). 

limit a,s n — > 00, 



For the bracelet graph Bn, 



the 



1. The threshold density C^c is the unique positive root 
o/C = I - ^e-^C (numerically, Cc = 2.496608/. 

2. The stationary density is 5/2. 

3. The final density /9„(A), as a function of initial den- 
sity X, converges pointwise to a limit p{X), where 



p{X) = min A 



-2A 



A<Cc 
A>Cc 





« 




FIG. 1: The graphs on which we compare and (c- the 
grid (upper left), bracelet graph (upper right), flower graph 
(2""^ row left), complete graph (2°"^ row right), Cayley trees 
(Bethe lattices) of degree d = 3,4,5 (S"^"^ row), and ladder 
graph (bottom). 



Part 3 of this theorem shows that despite the inequality 
Cs 7^ Cci a connection remains between the driven dissi- 
pative dynamics used to define Cs and the conservative 
dynamics used to define Cc' since the derivative p'(A) 
is discontinuous at A = <^ci the driven sandpile under- 
goes a second-order phase transition at density C.c- For 
A < Cc, the driven sandpile loses very few particles to 
the sink, and the final density equals the initial density 
A; for A > Cc it loses a macroscopic proportion of par- 
ticles to the sink, so the final density is strictly smaller 
than A. As Figure shows, the sandpile continues to 
evolve as A increases beyond Cc; in particular, its density 
keeps increasing. 

We are also able to prove that a similar phase transi- 
tion occurs on the flower graph, shown in Figure [T] In- 
terestingly, the final density p(A) for the flower graph is 
a decreasing function of A > Cc (Figure [2] bottom). 

Our proofs make use of local toppling invariants on 
these graphs. On the bracelet graph, since particles al- 
ways travel in pairs, the parity of the number of particles 
on a single vertex is conserved. On the flower graph, the 
difference modulo 3 of the number of particles on the two 
vertices in a single "petal" is conserved. 

One might guess that the failure of the density con- 
jecture is due only to the existence of local toppling in- 
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p C. =5/2 



2.5 



A 

->■ 2.48 



1 p 


: C=5/2 










A 

^^ ^ >■ 



Cc2.5 




= 5/3 



Cc 1-7 



FIG. 2: Density p(A) of the final stable configuration as a 
fiinction of initial density A on the bracelet graph (top row) 
and flower graph (bottom row) as the graph size tends to 
infinity. A phase transition occurs at A = ("c. At first glance 
(left panels) it appears that the driven sandpile reaches its 
stationary density C,s at this point, but closer inspection (right 
panels) reveals that the final density p(A) continues to change 
as A increases beyond C,c- These graphs are exact. 



variants, or else to boundary effects from the sinks. The 
ladder graph (Figured]) has no local toppling invariants; 
moreover, it is essentially one-dimensional, so the bulk of 
the graph is well insulated from the sinks at the bound- 
ary. Nevertheless, we find [s^ a small but undeniable 
difference between C,s and C,c on the ladder graph. 



CONCLUSIONS 

The conclusion of [5] that "FES are shown to exhibit 
an absorbing state transition with critical properties co- 
inciding with those of the corresponding sandpile model" 
should be re-evaluated. 

In response to this article, several researchers have sug- 
gested to us that perhaps the density conjecture holds for 
stochastic sandpiles even if not for deterministic ones. 
This hypothesis deserves some scrutiny. 

For the driven dissipative sandpile, there is a transition 
point at the threshold density of the FES, beyond which 
a macroscopic amount of sand begins to dissipate. The 
continued evolution of the sandpile beyond Cc shows that 
driven sandpiles have (at least) a one-parameter family 
of distinct critical states. While the stationary state has 
rightly been the object of intense study, we suggest that 
these additional critical states deserve further attention. 
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